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ABSTRACT 


The method of Time Domain Structural Synthesis is reviewed and examples of 
linear structural modification and sub-structure coupling are presented. The general 
formulation for both linear and nonlinear syntheses are compared to highlight the 
similarities between the two and the difficulties experienced with the introduction of 
nonlinear elements. Examples of linear and nonlinear sub-structural coupling are 
presented to emphasize the similarity of the procedures used. Finally, a general approach 
for the optimal design of nonlinear shock isolation for large structural systems is 
developed. The repeated solution of arbitrarily large finite element models with localized 
nonlinear components is made possible through the use of a highly efficient nonlinear 
transient re-analysis procedure. The method has demonstrated order of magnitude 
decreases in compute time over traditional direct integration methods and is exact, 
regardless of the nature of the isolation. An example is included which demonstrates the 


procedure. 
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I. INTRODUCTION 


With the introduction of Commercial-Off-the-Shelf (COTS) technology onto 
military platforms, comes the need to ensure that such equipment is properly isolated 
from shock. Since COTS are not designed to meet military specifications with regard to 
shock survivability, shock isolation systems must be developed which ensure equipment 
survivability in a combat environment. However, to physically test each equipment 
configuration until a suitable mount is found would become prohibitively expensive. A 
more cost-effective approach involves modeling equipment and mounting systems on a 
computer. Computer simulations could then be conducted to determine the feasibility of 
different mounting configurations. 

Using finite element (FE) techniques and powerful computers, it is now possible 
to construct mathematical models and conduct response analysis of complex structures. 
By "breaking down" a complex system into its most basic components (stiffness’, masses, 
and dampers), the behavior of the system may be determined through the solution of a 
series of less complex equations. Assembling the "pieces" of the model into a set of 
system matrices allows for examination of structural responses (static and dynamic) of 
the complete system. The final phase of this process, solution of the model responses, is 
the most computationally demanding step in the FE process. | 

Having completed a response analysis, engineers may wish to make modifications 
(small or large) to the system, thus requiring further simulation. Traditionally, this would 
require the engineer to repeat the entire process for each change desired. However, an 
optimization routine that would determine the optimal design for a given structure as well 
as perform the Finite Element Analysis (FEA) for each "design" would be far more 
preferable. Even this method is not the most desirable solution, though, since the entire 
model must be reconstructed for each change. 

A more efficient method, that requires a system model to be constructed only 
once, and allows for solution of relatively few system elements, has been developed. The 


method dealt with in this thesis involves the use of synthesis procedures to generate post- 





modification system responses. This synthesis method (Time-Domain Structural 
Synthesis (References [1] and [2]) uses pre-modification mode] data, structural 
modifications made (in the form of stiffness, mass or damping changes), and a series of 
equations to solve for the modified system response. In using this method, the engineer 
need only solve for the large, pre-modification model a single time. Once this solution is 
obtained, only the modification matrices need be developed. Specifically, time-domain 
synthesis requires the pre-modification system’s impulse response functions (IRF’s), and 
the convolution integral to solve for post-modification response. The governing equation 
of time-domain synthesis is a non-standard, non-homogeneous Volterra 
Integrodifferential Equation (VIDE) of the second kind. This VIDE is solved 
numerically in order to calculated the post-modification transient response. Thus, this 
technique allows for considerable reduction in processor time and improved capability 
for the application of optimization routines. 

Thus, while the FEA method requires a reconstruction of the entire model for 
each modification made, the Time-Domain Structural Synthesis method allows for a 
single, pre-modification model to be built. The response characteristics of this pre- 
modification model are determined one time. This response characteristic is then used 
repeatedly to solve for the transient response of as many mounting configurations as 
desired. The timesavings achieved over FEA method in this way alone may be 
considerable. Combine this with the fact that the FEA method requires solution of the 
transient response of the entire model while Time-Domain Synthesis requires solution of 
an arbitrarily small number of elements, and the foundation of optimization of large 


structures is built. 








ll. INTEGRAL EQUATION FORMULATION FOR TRANSIENT STRUCTURAL 
SYNTHESIS | 


The method of structural synthesis allows the user to calculate the dynamic 
response of two or more substructures coupled together. It also allows for solutions to 
systems in which individual structural members may have been added, removed, or 
modified in some way. As stated above, the motivation to use this method lies in the 
computational savings of solving only for the synthesized model rather than the, often 
much larger, completely assembled model. A time domain, structural synthesis allows 
for the reduction of a model such that only the coordinates of interest are retained.. It also 
allows for an exact synthesis of system damping. Perhaps most importantly, time domain 
synthesis allows for solution of the modified system using pre-modification transient 
response data. Thus, solutions of the modified system equations need not be obtained. 
Regardless of the modifications made to the original system, the same baseline, pre- 
modification transient responses are used to generate the synthesized response. This 
allows for more rapid solution of the multiple iterations often required in an optimization | 
routine. The material presented in this section is simply an overview of the background 
work of Reference [1] and serves as a foundation for objective of this thesis. This 
material is presented to give the reader a general understanding of the techniques used in 


synthesis routines, both linear and nonlinear, which follow. 


A. STRUCTURAL MODIFICATION 


Consider an arbitrary structure, as in Figure 2-1, which may be partitioned into 
two subsets of physical coordinates. For simplicity, examples programs and calculations 
were carried out using one and four degree-of-freedom (translation only), spring-mass 
models. Modifications were made only to spring stiffness (linear modifications only), 
with no changes in mass or addition of dampers. In each of these models those 
coordinates where modifications are installed are designated as Connection Coordinates 


and are represented by {x,}. Coordinates that are not directly involved in the structural 


5 








modification, but whose post-modification transient response are of interest to the user, 


are designated as Internal Coordinates, and are represented by {x;}. 


Figure 2-1: General NDOF System 


The following is an overview of the derivation of the equations used to perform 
post-modification response analysis (linear and nonlinear) after structural modification. 
A more detailed derivation is available in Ref. [1], and shall not be included here. 

Starting with the convolution integral, the total dynamic response of the system 


can be expressed as follows: 
X, (t) x, (t) e(/h.(t-t) h,,(t-—T) |/ EQ 
hae 7 ie ot . | i (t-T) h,(t- a iret oY 


Where the [h(t)] are matrices of impulse response functions (IRF’s), while the 
subscript h denotes the homogenous solution. The force vector {F} contains forces 
experienced by both the internal and connection DOF. Internal DOF will experience 
externally applied forces (denoted by superscript e) only. Connection DOF, however, 


will experience externally applied forces as well as reaction forces caused by the 
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modifications (denoted by superscript *) associated with the respective DOF. Or, in 


vector notation: 


oe F(t) 0 - 
F(p) IR@! LE Coe 


Thus, we can express the synthesized response as 


x; (t) - x; (t) h,.(t- 
pee “eof i (t- off (7 1) jdt (2.3) 


In the case of Equation (2.3), the superscript * denotes the synthesized quantity 
(displacement in this case). Also, note that the first term on the right hand side of 


Equation (2.3) is no longer simply the homogeneous response, but rather 
x, (t) x, (t) rj h.(t-t) h,(t-—T) ||/Fct 
ed ee +] — Ola (2.4) 
X(t) X.(t)), ol ha(t—7) h,.(t— 7) || F(t) 


which is a nonstandard, non-homogeneous Volterra Integrodifferential Equation 
(VIDE) of the second kind. . 

Once the synthesized, transient solution of the connection DOF has been 
determined, the transient response of the internal DOF may be calculated using the 


remaining portions of Equation (2.3). 





B. SUBSTRUCTURE COUPLING 


Still using Equation (2.1), where the subscripts c and i retain their original 
meaning, a governing equation for the solution of synthesized transient response under 
the conditions of substructure coupling shall be developed. 

Given two substructures, a and b (as in Figure (2.2)) the coordinates for each can 


be partitioned as 


{x (tyb=[x2 xt]! {x.(tybe[x? x?’ (2.5) 





Figure 2-2: General Substructure Coupling 





Just as with structural modification, internal DOF experience external forces only, 
while connection DOF experience the additional, reaction forces due to coupling. The 


force vector now becomes 


{FO }=[RHE (wf (2.6) 


where [R] (which reflects force equilibrium) is the Boolean matrix (each column 
of [R] is all zero except for a 1 and a-1 in rows corresponding to coupled coordinates). 


The synthesized response now becomes, after pre-multiplying by [R]', 


{x2} {z,(o}+ J [b.. ¢- OE © fart (2.7a) 
0 

where | 

{x (ty }=[RT {x,t 7 | (2.7b) 

and | 

[,.(t)|={RT"[h,. PR] (2.7¢) 


Finally, what is left, is a homogeneous, linear Volterra Integral Equation of the 


first kind for the coupling force vector. In order to generate a solution for {F (t)} two 


derivatives are taken, resulting in 


/h.. olf w}=-KF«}-|[h.c-ofE ohn (2.8) 





which may be solved numerically for the coupling force vector. This is 
substituted into Equation (2-3), that provides a solution for the synthesized, transient 


response of the coupled system. 


C. SOLUTION OF THE VOLTERRA INTEGRAL EQUATION 


A more detailed derivation of the solutions to Volterra integral equations is 


presented in References [1] and [3], a general summary is presented here. 


b(x). . 


h(x)u(x) =f(x)+ [K(x,6)u(Ede (2.9) 


Equation (2.9) represents the most general, linear integral equation in u(x). This 
becomes the Volterra Integral Equation when b(x)=x. Should h(x)=0 (not to be confused 
with the h(t) for IRF’s), this equation becomes the Volterra Integral Equation of the first 
kind. However, if h(x)=1, the result is a Volterra Integral Equation of the second kind. 

Since, as shown in Reference [3], any general initial value problem may be 
converted into a Volterra Integral Equation of the second kind, the solution of the 


response of a spring-mass-damper system readily lends itself to the Volterra Equations: 


d? d bs 
(IM)? + (CDE + CKD xc = F(t) (2.10) 


In order to solve a VIDE of the second kind, the integration interval is first 
divided into a set of n-subintervals with time step At=(x,-a)/n, where x, is the end point 
chosen for x. Setting to=a and t==at+jAt=to+At, the value of the solution at t; or x; will be 
referred toas u(t,)=u(x,)=u,;, f(x,)=f,. Similarly, the value of the Kernel K(x,t) 
at (x;,t)) will be K(x,,t,) = K, . Using the trapezoidal rule with n-subintervals, the 


integral is approximated as 











2 


5 K(x,t, )u(t,) + K(x, t, u(t, )+... 


[K(x,pu(tdt = At ; 
A + K(x, tut) +5 Ot, uct, ) 








QZ AN) 
t; —-a x-@a 
At = = , t Sx, jal x=x, =, 
j n 
which may be approximated by the summation 
* K( ty )uct )+K(x,t,) (t,)+ 
—K(x,t,)u ,t,)u bes 
u(x) =f(x)+At}2 07° °° se 
+...+K(x,t,_, )u(t,_,)+K(x,t, u(t, ) 
(2.12) 


Substituting for each time interval, simplifying, and transferring solution (uj) 
terms to the left and non-homogeneous (f;) terms on the right, the following triangular 


system of equations results. 


Uy Pe ‘es es a: =f, 

At ; At 

—Z Kioto + 1- Ky, u, ia ie. us =f, 
At At 

—F Bato — AtK.,,u, + i x, ho bes es = f3 
At At 

-Z Kao = AtK,,u, - AtK,.u, + a | =, 
At At 

—Z Botte — AtK ,u, - AtK ,,u, — es (1-8, Jo =f, 

(2.13) 





Note that the Kj, elements equal zero in Equation (2.13) due to h;(0)=0. Also note 
that the solutions are obtained by forward substitution, starting with ug=fp in the first 
equation being substituted into the second equation to find u;. The results of the second 
equation are used to find the solution of the third and so on until u, is solved. 


In the case of structural modification the governing equation is 


iM} co} 
fii0}=f00}- [bat i "Ife coy 


where, for a simple stiffness modification, the equation reduces to 


{x*(t) b= dx (phe ils. (t—1) KK" Hx" a Hac. (2.15) 


The {x,(t)} is the pre-modification system response, the IRF matrix [hi-(t)] 
represents the Kernel matrix, and the {x, (t)} vector is the desired, post-modification 
system response. As with the derivation above, the solution is obtained through forward 


substitution as follows 


{x*(t)}= [KERNEL] {x,(t)}. | | (2.16) 


As mentioned previously, structural modification already presents itself in the 
form of a VIDE of the second kind. Structural coupling, however, is in the form of a 
VIDE of the first kind. In order to use the methods described above, the governing 
equations for structural coupling must be reduced to a VIDE of the second kind. This can 
be accomplished by following the procedure outlined in Reference [3]. 

Thus, the Volterra Integral Equations of the second kind allow for a less 
complicated solution of nwieldy and complex matrix equations. The VIDE of the second 


kind becomes the basis of the Time-Domain synthesis solution. 
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D. EXAMPLES 


Example 1. One Degree of Freedom Stiffness Modification 


; k Ak 
\\\\ + —\W\\— 
y(t) 





x(t) 


| k+Ak 
=... AW 
| a 
y(t) 


x(t) 


Figure 2.3 One Degree Of Freedom Spring-Mass System 


Starting with the simplest possible modification, consider the single DOF spring- 
mass system shown in Figure 2.3. A spring-mass system is modified by increasing the 
stiffness (linear modification) of the spring to ground by some value, Ak. The 
synthesized transient response is generated by using the transient response of the original , 
system in conjunction with its impulse response function (IRF). Given a harmonic 
forcing function of the form FsinQt (where F is the amplitude and Q is the forcing 
frequency), the solution to this system’s equation of motion (mx + kx = Fsin Qt) is 


readily calculated as 
F QQ. 
x(t) = era Or - 2 Sno, (2.17) 


n 


and the general form for the IRF is 


11 

















2 OPO? 
h(t) = dYororet >; G95 yt sin(«,,t) (2.18) 


p=l p=r+l dp 


where: r - number of rigid body modes 
n - number of modes retained 
;°, 0; - the p™ modal response of the i", j® mode, respectively 
Wap -- the p'" damped natural frequency 
C, - the p™ damping factor 
Wap - the p'' natural frequency 


t - time in seconds 


which reduces to 





ieee (2.19) 
as 
The governing equation of the synthesis becomes 
Ak t 
x*(t) = x(t)- —|x"(t) sin(@, (t—1))dt. (2.20) 
n 0 


Since this is a convolution integral, the transient response will be calculated at 
times t=iAt, i=0,1,2,3,....,n. Letting Lj=sinw(GAt-jAt), the integral in Equation (2.20) 


can be replaced by a trapezoidal rule 


‘ Ak | 1 ‘ — — l i: ; 
x(t) = x(t) - At] 5 Lig (OAt) + >) Lx (jAt)+>Liox'GAt)} i =0,1,2,..,n 


}Fl 


(2.21) 


resulting in the lower triangular matrix as described in the previous section. The solution 


is obtained through forward substitution. Figure 2-4 provides a plot of the synthesized 
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transient response of the modified system and a comparison to the exact solution using a 
trapezoidal integration scheme. Note that, with respect to the resolution of the plot, the 
synthesized response exactly matches the exact response. A copy of the computer routine 


used is provided in Appendix A. 


Comparison of Synthesized and Exact Displacement VS Time 
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Figure 2-4: Comparison of Synthesized and Exact Displacement Response for 


Modified 1-DOF System 
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Example 2: Multi-Degree-of-Freedom Stiffness Modification 
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Figure 2.5 Multi-DOF System with a Stiffness Modification at DOF #2 


A somewhat more complex modification involves a multi-DOF system in which 
the modification (linear stiffness) and the DOF of interest are not connected, as shown in 
Figure 2.5. The spring-mass system is modified by increasing the stiffness of the second 
DOF by some value, Ak. With the second DOF modified, the transient response of the 
fourth DOF will be observed. Again, the synthesized transient response is generated by 
using the transient response of the original system in conjunction with its impulse 
response function (IRF). Given a harmonic forcing function of the form F(sinQt+cosQt) 
(where F is the amplitude and Q is the forcing frequency), the solution to this system’s 
equation of motion requires the solution of four simultaneous equations. Since the 
system transient response is not so simply obtained as in Example 1, a computer is used 
to solve for the pre-modification transient response. 


Applying Equation (2.18) once more, note that IRF calculation is simplified to : 


h(t) = 5 ES ete sin(w,,t) (2.28) 


=] dp 
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As in Example 1, the result is a convolution integral, which can be solved ina 
similar manner. Using a trapezoidal rule and forward substitution, the results are 
presented in Figure 2.6. Note that, as in Example 1, the synthesized solution is exact to 


within the resolution of the plot. 


Synthesized Response of DOF #4 
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Figure 2-6: Synthesized Transient Response of DOF #4 to a Stiffness 
Modification at DOF #2 
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II. NONLINEAR SYNTHESIS 


The procedure outlined below was presented previously (Ref. [2]) for conducting 
a linear synthesis. The background and formulation for the nonlinear synthesis process 
are similar and are taken from References [2] and [4]. Rather than restating the 
formulation presented in Reference [2], a quick overview will be presented. Examples of 
both linear and nonlinear synthesis shall also be presented to show that the synthesis 
technique is similar for each case (albeit somewhat more computationally demanding in 
the nonlinear case). This synthesis formulation combines the computational efficiency of 
the time domain structural synthesis with the ability to install nonlinear modifications to 
the pre-modification structure. Through an iterative process, the exact solution is 


obtained quite rapidly. 


A. GENERAL (LINEAR) SYNTHESIS FORMULATION 


In this introductory section, all discussion of synthesis formulation will focus on 
the less complicated, linear problem. An iterative solution procedure is developed for the 
linear problem, which readily lends itself to use in the nonlinear problem. Thus, by | 
understanding the formulation for the less complex, linear case, the reader is better able 
to understand the nonlinear problem. The section which follows deals specifically with 
the nonlinear synthesis formulation. 

Consider a general NDOF system, as in Figure 3-1, with at least one rigid body 
mode, to which structural modifications are to be made. The first step in this process 
involves the assembly of the system stiffness, mass and damping matrices, either 
manually or through some finite element modeling subroutine. When the system 
matrices are assembled, the system natural and damped frequencies, as well as the system 
mass normalized modal matrix ([®]) are calculated. With this data and a given span of 


time over which the system’s transient response is desired, the impulse response function 
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(IRF) matrix [h(t-t)] is determined for the pre-modification structure using the general 


IRF formula 


n 


Py? 
n(t)= So?ort+ >, uh sin(«,,t) (3.1) 
p=l 


p=r+l O gp 


where all variables are defined as for Equation (2-18). 


Figure 3-1: General NDOF System 


Now that the system’s pre-modification response function has been calculated, the 
system modification is "installed". Using the convolution integral equation (full 


synthesis equation), the modified system total dynamic response can be written as: 


t 


{x(t ¥ = {xp} + | [ha—a Roya (3.2) 


0 


where the force vector is now the modified force vector, which includes coupling 


forces between connection and base DOFs. Thus, the forces may be rewritten as: 
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F _— ES 
F= F™ hy = F™ — (AM..x 
R= ER + F = isa 7 [(Ac, (x, —y)+ AK, es = y)| 


+AC,x. + AK,x.) (3-3) 


c* 


where: 


F*-the coupling forces due to structure modification 


[AM], [AC], [AK]-the modification matrices 


* 


x",x", x *—— generalized (synthesized) responses after modification 


Combining equations (3-2) and (3-3) results in a nonstandard, nonhomogeneous 
VIDE of the second kind. The solution is obtained through the use of the trapezoidal 
quadrature rule as discussed in Section II, above. Also, by assuming no external 


excitations (i.e. {F.'}=0), the problem is reduced to the form 


X, i My 7 Fix. ; Frc. P Fax. 4 
X, = cc . 0 Fc, Fy ( * ) 


where [A] is defined as follows: 
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In this quadrature matrix, the subscript n denotes the number of time steps. The 


coupling forces of equation (3-4) are defined as: 


ieee oh 3-6.a) 
les [ verde ' _om 


fect [A pares} 2 


The solution to this problem is best obtained through an iterative procedure due to 


the large size of [A]. The iteration process followed is outlined below: 


1) Assume the force vector {F,} 


2) Calculate the response vector {xe}" 
3) Use {xe} to calculate {XV {k}, and a new {F;} 
4) Use the new {F_} and calculate a new {x.} 


5) Repeat steps 3 & 4 until fen ee - {xe} old < convergence tolerance 


Through the use of the trapezoidal quadrature and the iterative methods, 
approximations have now been introduced to the time domain synthesis method (which is 
an exact formulation in itself). Convergence and stability become issues of concern with 
the iterative procedure. Smaller time step size (At) will lead to an increase in stability 
and a more rapid convergence by keeping the norm of [A] small. With smaller time step 
size comes the need for larger computer memory requirements for a given time interval 
of interest. There may be limitations on the span of time which may be "observed" by 


this analysis, dependent on the computer system used, due to memory capabilities. 
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B. NONLINEAR SYNTHESIS FORMULATION 


Given that an entire model is rendered nonlinear by the inclusion of a single 
nonlinear element, traditional direct integration solutions become computationally 
expensive. By isolating the nonlinearities from the bulk model the majority of the model 
gains the benefits of linearity, simplifying calculations. Also, by retaining only those 
degrees-of-freedom (DOF) of interest to the analysis, a more efficient response analysis 
is developed. This formulation allows for repeated, rapid nonlinear transient analysis for 
large, linear structural FE models, which can have localized, generally nonlinear 
components. | | 

Using only those DOF of interest in physical (non-modal) coordinates, a transient 
analysis of an arbitrary structure may be conducted which is independent of model size. 
The retained DOF must include, as a minimum, those DOF associated with the nonlinear 
elements, and those DOF of which the synthesized transient response is desired. Thus, it 
is possible to solve for the transient response of an arbitrarily large model using an 
arbitrarily small number of DOF (limited only by the number of nonlinearities and points 
of interest). This feature has been demonstrated for both the linear and nonlinear 
formulation References [2] and [4], respectively. 

As mentioned previously, it is desirable to isolate the nonlinear elements from the 
bulk, linear model. In order to achieve this, the structure must be divided into a number 
of substructures, Figure 3-2. By dividing the structure in this way, linear portions of the 
model may remain linear. For those linear substructures subjected to excitation, the 
baseline transient response is calculated. Any substructure not subjected to this 
excitation only requires that the impulse response functions be generated. This is known 
as a substructure coupling approach. With the nonlinear elements isolated, this method 
provides computational advantage by exploiting inherent physical boundaries in the 
problem and maintaining substructure linearity. The synthesis is used to connect the 
linear substructures through the nonlinear elements, and to calculate the combined system 
nonlinear transient response. Once the nonlinear elements are connected, the entire 
system is rendered nonlinear and the nonlinear transient response is described exactly by 


the governing equations of synthesis. The iterative method introduced in the previous 
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section on general synthesis is applied to the nonlinear synthesis problem with similar 
results. One important modification, which simplifies calculations and reduces computer 
memory requirements, has been made in the solution process. Rather than generating the 
quadrature matrix at each iteration, the appropriate IRF for the point of interest is 
convolved with the associated force vector {F,} using the MATLAB CONV command 
(Reference [5]). For large systems, this results in a considerable reduction in the number 
of calculations required to complete the synthesis (and thus, improved timesavings). This 
savings in time and memory requirements allows for increased stability and convergence 
capability by permitting even smaller time step sizes for the same memory requirements 


as the general synthesis procedure described in the previous section. 





ft ye 


Figure 3-2: Substructure Coupling/Modification for Nonlinear 
| Synthesis 
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C. EXAMPLES 


The validity/accuracy of the time domain synthesis technique has already been 
shown in Refs. [3 and 4]. The following examples are provided simply to demonstrate 
the ability of this method to solve for large, linear and locally nonlinear models with the 


same algorithm and equivalent performance. 


Examples 1 and 2: Mass-Plate System 


Consider the following mass-plate system, presented as Figure 3-3: 


ea 
You} I k, 





Figure 3-3: Mass-Plate System Experiencing Base Excitation 


This mass-plate system is the same as used in Ref. [4]. The stiffeners and 
dampers are added as modifications to the baseline system. The same system is used in 
both Examples 1 and 2. However, in Example 1, the stiffeners added possess a linear 
force-displacement characteristic of 500 lb;/in, while the stiffeners in Example 2 possess 
the nonlinear force-displacement characteristic curve shown in Figure 3-4. In both 
examples, the dampers possess a linear damping characteristic of 2 lby-s/in. The 
stiffeners used in Example 2 are representative, but not an exact reproduction, of the 
starting point to be used for the optimization routine. 
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The original structure is represented by the plate, which has spring-to-base 


attachments at each corner, and the attached computer which is located off-center, above, 
and connect to the plate with a linear (1000 Ib;/in) spring. (Note that base isolators are 
not installed in the pre-synthesis model, making it a free-free structure.) The plate- 
computer system is approximately 51,500 DOF. The computer 1s modeled as a single 
lumped mass with a single DOF in translation only. The base excitation, y(t), is a blast 
function and is presented as Figure 3-5. 

Figure 3-6 is a plot of the transient response of the computer to a blast function 
with linear isolators installed. Figure 3-7 is a plot of the transient response of the | 
computer using the same base excitation with nonlinear stiffeners and linear dampers 
installed. The synthesis algorithm for both analyses is the same and is provided in 
Appendix B. Note that the only change made to the algorithm involves replacing the 
linear isolator with the nonlinear isolator. The synthesis performed for the linear isolator 


required 22.09 sec, while the synthesis for the nonlinear isolator required 50.84 sec. 


Nonlinear Spring Force VS Displacement 
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Figure 3-4: Nonlinear Spring Force-Displacement Characteristic 
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Base Excitation (Blast Function) VS Time 
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Figure 3-5: Base Excitation Y(t) (Blast Function) 


Equipment Response VS Time (Linear Spring) 
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Figure 3-6: Transient Response of Computer with Linear Spring and Damper 


Modification 
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Equipment Response VS Time (Nonlinear Spring) 


Displacement (in) 





0.02 
Time (sec) 


Figure 3-7: Transient Response of Computer with Nonlinear Spring and 


Linear Damper Modification 
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IV. OPTIMIZATION 


The background information presented here is an overview of what is presented in 
reference [6]. With the introduction of Commercial-Off-the-Shelf (COTS) technology 
into military platforms comes the need to ensure that such equipment is properly isolated 
from shock and vibration. Since COTS are not designed to meet military specifications 
with regard to shock survivability, shock isolation systems must be developed which 
ensures equipment survivability in a combat environment. However, to physically test 
each equipment configuration until a suitable mount is found would become prohibitively 
expensive. A more cost-effective approach involves modeling equipment and mounting 
systems on a computer. Computer simulations could then be conducted to determine the 
feasibility of different mounting configurations. | 

The problem is still not so simple as it may appear. Design of a shock only or 
vibration only - isolation system may be relatively easy. However, when trying to isolate 
both shock and vibration, the optimal design for one is in conflict with the other 
(References [7] and [8]). This seriously complicates the optimal design of a shock and 
vibration isolation system. Thus, the objective of the optimization becomes the 
minimization of some aspect of the system response while also constraining other 
response characteristics within some reasonable bounds. 

Even this does not fully describe the complexity of the problem presented. Since 
the system to be isolated is typically large, the calculations alone are demanding. Add to 
this the need to explore variations in a large number of design variables and the 
optimization becomes more cumbersome. Finally, given that an entire model is rendered 
nonlinear by the inclusion of a single nonlinear element, traditional direct integration 
solutions of a single configuration become computationally expensive. Combine this 
with an optimization routine requiring numerous configuration tests, and the optimization 
of nonlinear isolators becomes, for all practical purposes, prohibitively expensive. 


However, by further developing the work presented in Refs. [3 and 4], an efficient 
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process for determining the optimal nonlinear isolation characteristics for large, locally 
nonlinear system responses can be obtained. 

By isolating the nonlinearities from the bulk model the majority of the model 
gains the benefits of linearity, simplifying calculations. Also, by retaining only those 
degrees-of-freedom (DOF) of interest to the analysis, a more efficient response analysis 
and optimization routine is developed. This formulation allows for repeated, rapid 
nonlinear transient analysis for large, linear structural FE models, which can have 


localized, generally nonlinear components. 


A. OPTIMIZATION ROUTINE FORMULATION 


In the formulation of an optimization routine, several specific tasks must be 
accomplished to define the optimization space. First, an objective function must be 
defined which specifies the quantity to be minimized or maximized. Second, design 
variables (those quantities to be varied in the optimization) must be determined and 
initialized. Finally, the problem must be constrained in some way, and constraint 
functions, which are themselves functions of the design variables, must be developed. 
These constraints help limit the size of the optimization space, allow for the imposition of 
physical restrictions, and when properly chosen, can reduce the optimization search time. 
Three types of constraints may exist in an optimization problem. First, an inequality 
constraint is one in which the value calculated from the constraint function is less than or 
equal to some prescribed limit. Second, an equality constraint is one in which the 
constraint function value is equal to some prescribed value. Finally, a side constraint is 
one where some upper and lower limits are placed on the design variables. 

For the design of an isolation system, the minimization of the system’s response 
due to some excitation may be the ultimate objective. The minimization of response 
history away from equilibrium, the absolute or relative displacement, or perhaps the 
absolute acceleration of the system may be important. Any one of these responses may 
become the objective function while the other responses may be used as constraints on 


the system. While minimizing the maximum acceleration experienced by the equipment 
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may be the primary concern, returning the system to equilibrium as quickly as possible, 
and the absolute displacement or the stretching of the isolators (relative displacement) 
must fall within some limits to avoid undesirable damage. To achieve the objective, 
physical parameters such as isolator stiffness and damping might be altered. Thus, as in 


Reference [2] the following optimization problem is generated: 


Minimize (Objective Function): 

Maximum system acceleration due to base excitation 
Subject to (Constraints): 

Maximum dynamic isolator stretch/compression < limit 

Motion away from equilibrium < limit 

Maximum equipment acceleration < limit 

Min and max changes in isolator stiffness < limit 

Min and max changes in isolator damping < limit 
Design Variables: 

Changes in the isolator stiffness characteristics 


Changes in the isolator damping characteristics 
1. Calculation of Objective Function 


The objective function is based on the system’s (computer) maximum acceleration 
due to some base excitation, which may be any random vibration or shock input. The 
displacement response is calculated first, then a finite difference scheme (which accounts 
for beginning (forward scheme) and ending (backward scheme) history points) can be 
used to calculate the acceleration response. From this acceleration time history, a 
maximum of the absolute accelerations experienced is readily obtained. By minimizing 


the accelerations experienced by the computer, the potential for damage is minimized. 


| (4-1) 
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The objective function becomes 


f= max 


d’x(t) 
dt? 

















For each iteration of the optimization process Equation 4-1 must be evaluated. 
This is accomplished by first calculating the system (computer) response through the use 
of the Time Domain Structural Synthesis for each variation of the optimization routine. 
The transient response of the computer is generated then used to evaluate the objective 


function. 


2. Calculation of Constraints 


The first two constraints noted are of particular concern in the formulation of this 
problem. Possibly the most important constraint is held on the maximum dynamic 
isolator stretch/compression. While the isolator stiffness (and therefore resistance) 
increases with increasing stretching or compression, there are physical limitations to the 
amount of stretching or compression an isolator may sustain without damaging or altering 
its force-displacement characteristic curve. Therefore, it is desirable to constrain the 
relative displacement between the plate and the base to some reasonable value within 
equipment specifications. In order to calculate this constraint, the transient response of 
both the base excitation and the structure (plate) at the points of isolation must be known 
(or calculated) for the time history of interest. Since the base excitation is prescribed, 
only the structure’s displacement time history need be determined. The time domain 
structural synthesis (nonlinear) technique is applied to calculate this displacement time 
history. 

The second constraint is a limit on the time required for the system to return to 
equilibrium after some base excitation. Since the primary concern, for this problem, is 
the isolation of equipment in a combat environment, the assumption is made that it is 
desirable to stabilize the equipment as quickly as possible without introducing large 
accelerations. This constraint, if overly restrictive, can conflict with the objective 
function. In order to return the system to equilibrium rapidly, larger accelerations will be 


developed than if the system is brought to equilibrium more gradually. Thus, there is a 
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trade-off between ensuring the equipment is protected from violent accelerations and 
returning the system to equilibrium rapidly. Thus, with the objective function taking 
priority, this constraint is simply a limitation on the amount of time the system requires to 
return to some acceptable level of stability. Note, that in many environments (seismic 
excitation of a structure, for example) a rapid return to equilibrium may not have any 
significance, and would not be considered a constraint. 

The third constraint helps drive the optimization routine past any potential local 
minima which would produce an undesirable level of acceleration. This constraint limits 
the maximum acceleration that may be experienced by the equipment and ensures that 
local minima which do not meet this criteria are discarded as invalid solutions. 

The last two constraints noted are simply limitations on the degree to which 
values of the coefficients of the characteristic force-displacement curve of the isolator 
may change for each optimization iteration. These constraints can be controlled by pre- 


defining design variable change limits for the optimization function. 
3. Design Variables 


Isolator stiffness and damping characteristics are changed during this procedure 
to optimize the response of the equipment. However, there are limitations to the amount 
of change these characteristic curves may undergo. The most obvious limitation is that 
the lower bound for the isolator stiffness characteristic is not equal to the zero. Use of 
this stiffness would be equivalent to the complete removal of the isolator. The upper 
bound of isolator stiffness and damping is governed simply by what is physically 
reasonable. One additional limitation is introduced with the nonlinear problem. As the 
characteristic force-displacement curve turns to change slope with increasing 
displacement/compression (Regions I and II in Figure (4-1) below), the new slope must 
be greater than or equal to zero. A negative slope in these regions will result in an 
instability in the synthesis procedure, which may results in an inability of the synthesis 
method to converge. A negative slope in these regions also indicates a system which is 


physically infeasible. In the case of a negative slope, spring forces would not oppose 
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compressions/expansions, but rather be sympathetic to them. This problem can be 
avoided by proper choice of the lower and upper bounds set on the design variables 
which govern slope in regions [ and III, such that a negative slope is not possible without 
violating boundary conditions. By properly choosing two points (or two slopes) on a 
force-displacement space and reflecting them across the displacement axis, a force- 
displacement characteristic curve is generated. Allowing the optimizer to alter the 
points/slopes chosen (such that the constraints prevent a negative slope in Regions I or Il 
as mentioned above), a force-displacement characteristic of the optimal calaiers can be 


generated. 


B. OPTIMIZATION COMPUTER CODE 


The optimization tool used to solve this problem, as in Ref. [3], is the MATLAB 
optimization toolbox 1.0d. Since the problem posed is nonlinear, constrained, and 
multivariable, the function CONSTR.m is used since it is specifically designed to solve 
such problems. CONSTR uses the Sequential Quadratic Programming (SQP) method in 
order to solve the problem (Reference [9]). The SQP method is used to determine the 


search directions for the design variables for each iteration (Refs. [5] & [9]). As with 


many nonlinear optimization computer routines, this method has the important limitation . 


that the optimal solution found might only be a Jocal solution. Thus, improvement might 
be made upon the solution obtained. This highlights the need for the user to have a good 
physical understanding of the problem so that an appropriate starting point is chosen. 
Another difficulty, which might be encountered, occurs when the problem posed is 
infeasible. In this case, CONSTR.m attempts to minimize the maximum constraint 
value. This can lead to the search pattern extending outside the defined optimization 
space, the user defined iteration limit being exceeded, and/or no optimal solution being 


found. 
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Figure 4-1: Sample Nonlinear Force-Displacement Curve 


The computer code used in this optimization routine follows the flow path 
presented in Figure 4-2. Model characteristics are loaded from a NASTRAN data file 
through the use of the function fReadNASModesPCH.m. Model natural frequencies 
{@n}, Impulse Response Functions (IRFs), and the forcing function are generated prior to 
starting the optimization routine. A starting point and upper and lower bounds for the 
optimization design variables are defined and CONSTR is invoked. For each iteration of 
CONSTR, the optimization function evaluates the post-modification model (based on 
design variables), calculates the value of the objective function, and verifies that 
constraints have not been violated. The routine remains within the optimizer and 
continues to modify the model until an optimal solution is located. A copy of the 


computer code used is enclosed as Appendix C. 
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Load NASTRAN Data 
(fReadNASModesPCH.m) 


Calculate IRF/ fn 
(f1RF.m ) 


Generate Forcing Function 


(fBlastForcing.m ) 


Define Initial Optim- 
ization Point 


Commence Optimization 
(CONSTR.m optfun.m) 





Calculate Post-modification 
Response (fSynth.m) 


Modify Design Variables 
as Reqd (optfun.m) 













Calculate Objective Funx 
and Constraints (optfun.m) 


Synthesis Convergence Test 
(fConvergeCheck.m) 





Optmal 
Solution 


Save and Plot 
, Results (fPlotter.m) 





Figure 4-2: Optimization Routine Flowpath 
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C. EXAMPLE 


Consider the following computer-plate system: 


¥ (DEY 0671-300) 
Youl re i, 
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Figure 4-3: Plate-Mass System Experiencing Base Excitation 


The original structure is represented by the plate, which has spring-to-base | 
sencheaents at each corer, and the attached computer which is located off-center, above, 
and connect to the plate with a linear (1000 Ib//in) spring. (Note that base isolators are 
not actually installed in the pre-synthesis model.) The plate-computer system is 
approximately 51,500 DOF. The computer is modeled as a single lumped mass with a 
single DOF in translation only. The base excitation, y(t), is a blast function identical to 
that presented previously in Figure 3-5 | All solutions are calculated using a MICRON 
266 MHz Millennium computer. 


35 





1. Problem Formulation 


The following is the formal problem statement for optimization: 


k, _stretch a Gee eqs 
51 = max k, _ stretch aan eae ai 


k ._ stretch 
~ maxk,_stretch 


ie 
2! 


ee te ee 
* 2 
nox{ fh (1) | 


d?x(t) 
dt? 





Minimize:  F(Ak,,Ak.,Ac,,Ac.) = mex 








g, ~1<0 7.21<x,, $9.65 


x"(0)| dt 





Subject to: 2, = 450 <k,, < 700 


w—Xen®COl cy 01S kK. <710 
ea Tmax x,_accel ee 
xX. dynam 
ge <0) 00<c, $5.0 


max x,_ dynam 


where: y,(t)’ = translational motion time history for the computer 


k,_ stretch = isolator stretch between base and plate 

k,_stretch = isolator stretch between plate and computer 

max k,_stretch = max. allowable isolator stretch between base and plate 

max k,_stretch = max. allowable isolator stretch between plate and computer 
X-_accel = absolute acceleration of computer due to shock [in/s”] 

max X,_accel = max. allowable absolute acceleration of computer [in/s”] 
X-_dynam = absolute dynamic displacement of computer due to shock 

“max x,_dynam = max. allowable dynamic displacement of computer 

Xb1 »Xb2 Kb1 >Ky2 = plate isolator design variables 


Cp = plate damper design variable 
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The following information is provided as a starting point for the optimization 
search routine. These points produce a rough approximation of the force-displacement 
characteristic curve presented in Figure 4-4. Table 4-1 provides these points along with 


upper and lower limits for more reasonable isolator profiles. 


Nonlinear Spring Force VS Displacement 
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Figure 4-4: Starting Point Force-Displacement Characteristic 


Kp Kp2 Cb1 


1.50 


Upper 10.00 700 710 5.0 


Table 4-1. Optimization Start Point Inputs 


C7 











The following inequality constraints are also placed on the optimization routine: 


max k, stretch < 0.50 inch 
max k,_ stretch < 0.50 inch 
max x,_ dynam < 1.00 inch 


max x,_accel < 15 ¢’s 
As mentioned above, the proper choice of the starting point is of particular 
importance. Since the optimization routine solution may result in a local minimum, the 
user must use his/her knowledge of the system to choose an appropriate starting point. 


2. Optimization Results 


This example is executed using the following main programs and their associated 


subroutines: 

Main Programs: Subroutines Called: 

1. f{MainOpt fReadNASModesPCH 
fIRF 
fBlastForcing 

2. Optfun fddot 
trapz 

3. Synth2opt fSymmetricStore 
conv 
fConvergeCheck 


fExceed_Iter 
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When the optimization routine is complete, the following information and plots 


are generated which help the user understand the optimization process. Table 4-2 is the 


display presented on the MATLAB Command Window, which provides information 


regarding the progress of the optimization routine. The first plot, Figure 4-5,isa 


comparison of the starting point acceleration response to the "optimal design" 


acceleration response of the computer after the modifications have been installed. The 


starting point design experiences a maximum acceleration of 23.4 g's while the optimal 


design experiences a maximum acceleration of 14.3 g’s. 


f-COUNT FUNCTION 


6 
12 
18 


24° 


30 
31 


23.9362 
14.7974 
14.7845 
14.3508 
14.3472 
14.3472 


MAX{g} 
7.93622 
-0.12069 
-0.120815 
-0.125029 
-0.125064 
-0.125064 


Optimization Converged Successfully 


ee oe ee 


Procedures 


Hessian modified 


Hessian modified ' 


Hessian modified 


Table 4-2: Optimization Routine MATLAB Generated Output 
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Acceleration VS Time (Comparison) 
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Figure 4-5: Comparison of Acceleration VS Time Between Starting Point and Optimal Design 


Modifications 


Figure 4-6 is a comparison plot-of the starting point displacement response to the 
optimal design displacement response of the computer. Note that while maximum 
displacements have not been altered to any great degree, there 1s a more gradual change 
of displacement with time tn the optimal design, leading to lower accelerations. 


Effectively, the natural frequencies of the system have been lowered. 
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Displacement VS Time (Comparison) 


0.8 X--Original Design 
__--Optimal Design 
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Figure 4-6: Comparison of Displacement VS Time Between Starting Point and 


Optimal Design Modifications 


Figure 4-7 is a plot of the slopes that characterize the nonlinear spring force- 


displacement profile at each iteration of the optimization routine. Note that the slope of 


regions J and III become negative in iterations 17 and 23. These points violate a 


constraint, are not valid solutions, and are thus rejected by the optimization routine. 
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Characteristic Stiffener Slope VS Iteration (Base Stiffener) 
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Figure 4-7: Change in Stiffener Slope (Base Isolator) VS Iteration Number 


Figure 4-8 is a plot of the installed damper characteristic slope versus iteration 


number. 
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Characteristic Damper Slope VS Iteration (Base Damper) 


Computer Damper Slope (#-s/in) 
(a) 
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Figure 4-8: Base Damper Characteristic Slope VS Iteration Number 


Figure 4-9 is a plot of the value of the objective function versus the number of 


iterations. 
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Objective Function Value VS Iteration Number 


a Stang ncn aw es Mocee donc veseeee peceebouserscce desccoee ncdeercckeecencbustecesleceeceteecesessbaccoucesecesws 


= 
@ 


seas eeeeneaweoue 80804 OEE Se EEE HME SSSA EERO EE EPR REE EHS RES Peewee eae 


Objective Function Value 
© 


ab 
~J 


eee Cees Pee eee ee eee ee ee eee ee ee eee eee ee ee eee ee eee eee ee eee ee eee eee ee eee 





0 5 10 15 20 25 30 35 
iteration Number 


Figure 4-9: Objective Function Value VS Iteration Number 


Figures 4-10 and 4-11 present the optimal design acceleration and displacement 


response time histories, respectively. 
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Figure 4-10: 


Computer Displacement (in) 


Figure 4-11: 





0 0.005 0.01 0.015 


“0 0.005 0.01 0.015 


Acceleration VS Time (Optimal Design) 


. 
Pe tell eel Pee ee) ee el 
’ * Py 


2 ee eee ene Ant 


eee ers eee) ee ee kd cia eee eee eee 












0.02 0.025 0.03 0.035 0.04 


Time (sec) 


Optimal Design Acceleration VS Time Response 


Optimal Design Displacement VS Time 


. 
ee eo cere wee eee OHH OOO Ee OE EEE DESC REESE SSSR DETER SRE ESSE SESS KeSenassaBewsrEsaaBenaTe es 
rf . . * 


wenenamaeauhostansaee eee em abe ORE OS EROS Ee SOOM SESE SESS TSE ESE SHS OTHE SEE SASE EHS SES 


. « * . 
ee re ea Oe eR ESSE ETCH TERE SEE HE Hee ES ORT Ema TeRewase see ssonaranehaseneanssanavesenhasen 
. 


Been cee mee tw coc ES ee Se NO i nt eae enn a = oom 
ewe ORO w eee ee eOee ene ee meee e eee aes hegeere sr Hesse Began seasereeEsweseZenesasasaneerEaaaane 





0.02 0.025 0.03 0.035 0.04 


Time (sec) 


Optimal Design Displacement VS Time Response 


45 








Table 4-2, and Figures 4-5 and 4-6 give a clear picture of how modifications to 
the structure can have such a large impact on system response. Since the objective 
function is to minimize the maximum of the absolute of the acceleration, it is readily 
visible, in Figure 4-5, how much the design has improved from starting point to "optimal" 
point. This reduction in the value of the objective function is even more apparent in 
Figure 4-9. Starting with an initial value of 23.4 g’s (with the starting point design 
variables), the optimization routine quickly reduces the value of the objective function by 
nearly 40%, to 14.3 g’s. 

Overall, this optimization required 31 iterations to determine optimal values for 
the design variables to minimize the objective function and meet the constraint 
conditions. This entire process required only 15 minutes and 47 seconds (for an average 
of 30.5 seconds per iteration). Possibly somewhat more impressive is the fact that an 
optimal solution was found so rapidly given that 5 design variables needed to be 
manipulated. Thus, this optimization routine calculated the transient response of a 
51,500 DOF structure, determined the maximum acceleration experienced for each 
modification and verified that constraints were met, given 5 variables to be manipulated, 
and still produced a valid solution in roughly 16 minutes. 

To verify the validity of the solution, several starting points close to the "optimal" 
were tested. Each time the same final solution was obtained in less than 10 minutes. 
Table 4-3 presents the starting points and time requirements for the optimization to return 
a solution. In each case, the optimal values for design variables are identical to those 


presented in Table 4-2 and are not reproduced here. 
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RUN # Starting Point Values Time Req’d 


a a 


3 1.35 9.55 490 704 4.9 375.90 


(sec) 


Table 4-3: Comparison of Starting Point Values VS Time Required 
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V. CONCLUSIONS/RECOMMENDATIONS 


This thesis was undertaken to investigate the feasibility of applying an 
optimization routine to the isolation (from shock and/or vibration) of large, locally 
nonlinear structures. Traditional FEA methods require the user to rebuild the model for 
each modification and the solution of the entire model. However, using the Time- 
Domain Structural Synthesis method, the pre-modification model is solved a single time 


and the solution of the transient response of an arbitrarily small number of DOF is 


required, resulting in exceptional timesavings. The work presented in References [1], [2], 


and [4] is the foundation for this thesis and has been presented, in part, as background. 
An important modification to the methods presented in these references, which improves 
time and memory requirements, was used in this thesis. While this modification did not 
significantly alter the procedures used, it removed the need to generate the quadrature 
matrix at each iteration and resulted in a more rapid synthesis solution. This 


improvement was incorporated in the synthesis technique and coupled with the 


optimization routine to determine the optimal configuration of an isolation system given a 


pre-modification model. 


A. CONCLUSIONS 


The following conclusions can be drawn from the analyses presented: 


1. The time-domain structural synthesis method is an efficient and accurate 
method of calculating the transient response of a structure following 
modification. The position and degree of modification is unrestricted by the 
method. Since numerical techniques are used to produce a solution, some 
error is introduced, through approximation, to the final solution. Errors 
introduced by numerical solution are limited by decreasing the time step size 


used. By generating only those values used for each convolution, at each 
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iteration, time step size may be further reduced since less memory is required 
(compared to generating the quadrature matrix for each iteration). 

Use of the time-domain synthesis method is integral in the optimal design of a 
nonlinear isolation system for a large structural system. Having demonstrated 
order of magnitude (or better, Reference [4]) improvements over traditional 
FE methods, in calculation time, the result is more rapid solution of the 
optimization problem. Through the use of this technique an optimal isolation 
system for a 51,500 DOF structure was generated in less than 16 minutes 

(15 minutes 47 seconds), using nonlinear stiffeners and linear dampers as 
design variables. After 31 different configuration tests (averaging 31 seconds 
per configuration) an optimal solution was presented. As a comparison, tests 
in Reference [4], using the same model, showed that a direct FE solution 
required 43 minutes 41 seconds to solve for the transient response of a single 
configuration. Further, the older synthesis routine used in Reference [4] (in 
which the quadrature matrix is generated for each iteration) required 4 
minutes 15 seconds to produce a transient response of a single configuration 
of the same model used in this thesis. Recognizing that the processor used in 
Reference [4] is different than the one used here, a better comparison of time 
requirements must be based on a comparison of the processor speeds used in 
each test. Based on the different processor speeds used, there is still better 
than an order of magnitude reduction in calculation time over the direct FE 
method, and an 80% reduction in calculation time over the older synthesis 


routine . 
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B. RECOMMENDATIONS 


While conducting this study, several important areas for improvement were noted. 


Recommendations to address these potential improvements include: 


e Develop a more intricate nonlinear stiffener generator in the optimization 
routine and a more robust decision algorithm for the synthesis to determine 
stiffener force/deflection as the synthesis converges. 

e Develop a more complex and realistic model of the damping portion of the 


isolator. 
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APPENDIX A. TIME-DOMAIN SYNTHESIS COMPUTER CODES 
FOR EXAMPLES 2-1 and 2-2 


The following is a brief description of the MATLAB computer codes used to 
perform the calculations necessary for Examples 2-1 and 2-2. Both examples shown 


present a linear structural modification. 


e VIBEIDOF.m - performs structural modification on a 1-DOF spring-mass 
system. Using time-domain structural synthesis to calculate the synthesized 
transient response of the post-modification system to a base excitation. 

e VIBEMDOF.m - performs structural modification on an M-DOF spring mass 
system experiencing base excitation. Allows the user to. choose the location 
of the modification and the DOF where the transient response is desired. This 
program also uses the time-domain synthesis technique to calculate the post- 


modification transient response. 
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VIDE1DOF .m 


KREKEKKKEKEKKERE KEK KEKE KK KEKE KEKE KEKEEK 


oe 


oe 


Define constants 
KKK KKK KKK KKK KKK KEKE KKK KEKE KEK 


fa) 


oe 


kspring=100.0; 
mass=100/386.4; 
dkspring=0.4*kspring; 


spring stiffness #/in 
mass (#m) 
spring modification 


oe 


dP 


W=250; % Forcing Frequency 
1 0 Fad a % Forcing Amplitude 
zeta=0.05; * Damping Factor 


KR KK KKK KKK KR KKK RK RK RK RK KK KEKE KE 


oe 


de 


Set time step information 
KKEKEKKEKKEKEKKEKEKEKEK KK KKK KE KKK ERE ESE 


de 


b1=0.0- $Start time 

tr=.10; Stop time 

dt=le-3; | S6Time step size 
numstep=(tf-ti)/dt; Number of time steps 


KEKKKKKKEKKEKEKKKEEKEKKEKKEKEKEKKKEKEKEKKEEKE 


de 


dO 


Calculate several constants 
KKK KEKE KKK KKK KEKE IERIE KEES 


oo 


Natural Freq. (rad/sec) 
Damped Nat Freq 
Exact Soln Consts 


de 


whn=sqrt (kspring/mass) ; 
wd=wn* sqrt (1-zeta%*2) ; 
cl=zeta*wn; 

e2=ci* 2. 
c3=c2+wd*24+2*wd*w+w’ 2 ; 
c4=c2+wd*2-2*wd*W+W%2 ; 
c5=W*2-c2-wd%2; 
c6=W*2+c2-wd“%*2; 


fe) 


oe 


KREKKEKKKKKEKKEKKKKKREKE KKK KKE KEKE EK KKK KEKE KEK EK KEKE KE KKREK EEE 


de 


oe 


Initialize Kernel Matrix and begin calculations 
RKEKKKEKEKKEREKKKEKRKEKEKEKREKE KEK KEKE KK RK RE KK KEKE KE KK KEKKERKEEERE 


de 


K=zeros (numstep+1,numstep+l1) ; 
Rly La 0e 
C=acE- 
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for 1=2:numstep+1 


ct2=(1-1) *=dt; 

£6 Gave 
tau=(j-1) *dt; 

if 2S]  S$Lower Triangular 
if j== 6First Column 


kernel=(1/wd) *exp (-zeta*wn* (t-tau) ) *sin(wd* (t-tau) ); 
K(i,j)=C/2.*kernel; 


else %2nd + columns 
kernel=(1/wd) *exp(-zeta*wn* (t-tau) ) *sin(wd* (t-tau)); 
K(i,j)=C*kernel; 


end . 
elseif 1==j SMain Diagonal 


kernel=(1/wd) *exp (-zeta*wn* (t-tau) )*sin(wd* (t-tau) ); 
K(i,j)=1.+C/2*kernel ; 
end 
end 
end 


KK KK KKK KKK IK KK KK KR KKK KER KER KKK KKK KKK KEE KEKE KEKE KEKE EESE 


Calculation and assembly of the x(t) vector 
ee KK KKK KKK KKK KR KKK KK EKER K KKK KKK KKK KKK KEKE KEKE 


dP oP dO 


Cp=c1 20 ts; 

tp=tp’; 
pl=f0.*(c5.*sin(tp.*W)+2.*c1.*W.*cos(tp.*W) )/(c3+c4) ; 
p2=f£0.*W. *exp(tp.*(-cl)) .*(c6.*sin(tp. *wd).. 


+2.*cl.*wd.*cos(tp.*wd)) ./ (wd*c3*c4) ; 
x=pl+p2; 


xstar=K\x; 
maxK=max (eig(K) ) 


wnmod=sqrt ((kspring+dkspring)/mass); %SMod nat freq 


KK KK KK KK KK KK RK KEKE KKK KEKE KKK KEKE KK KEKKKREKE KEKE 


Recalculated constants for new wn 
KKK RR RRR RE KEK KEKE KKK KKK RK K KKK KKK KKK RK KK KK KKK 
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dP dP dO 








clil=zeta*wnmod;' 
G22=cli*o; 
C33=c22+wd*2+2*wda*w+w2 ; 
C44=c22+wd*2-2*wd*Ww+w’ 2 ; 
c55=W*2-c22-wd“%*2; 
c66=W%24+c22-wd%2; 


pi1l=f£0.*(c55.*sin(tp.*W)+2.*c11.*W.*cos (tp. *W) )/(c334+c44); 

p22=f0.*W.*exp(tp.*(-cl11)).*(c66.*sin(tp.*wd)+2.*c11.*wd. *... 
cos (tp.*wd))./ (wd*c33*c44) ; 

xtrue=pl1l1+p22; 

dd +; 1:)=x- 

dd(:,2)=xstar; 

dd(s,3)=xtrue; 

£=0.05.*sin(tp.*tp); 


plot(tp,x,'-x',tp,xstar) ; 
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VIDEMDOF.m 


LEKKKKEAEKHKA KKK KKK KEKE ERK EK KEE KER ERE RE 


%¢ Enter Time Step Info and other constants 
LEKLERER REAR RARER RK AR RRR AIK RARER A AK SS 


ti=0..0; % Start time 

ab=000L; % Time step size 

ti=.50; % End time 

numstep=(tf-ti)/dt; % Number of time steps taken 
W=24.0; % Forcing Frequency in Rad/sec 
YO=s=i..0; % Forcing Amplitude 

nlin=1; % Activates NL Spring 


Sok kk ee A a ee He eR RK KK KK ERK KK KK KK KE 


% Define stiffness, mass and damping 
Sock kk RE aK KI RR KKK KK KK KK KK KK A KEK KK 


k=10 .250)- Stiffness vector 
m= (01. .02)3 Mass vector 
Z=([0: OQ) .Damping ratio vector 


K=[k(1)+k(2) -k(2);-k(2) 
M=(m(1) 0;0 m(2)]; 
dof=length(m) ; 


(2)]; % Define [K] 
Define [M] 
define # of system DOF’s 


dP doe AK’ dO oP AO 


Qk ke Ra eR RK KR RK KKK RR KKK KK HK KKK KKK KK 


dP 


calculate [PHI] and wn’s, wd’s 
Se KK KK KE KK KKK KKK KEK KEKE KEKE KKK KE KKEKEE KEKE KEE 


ae 


[omega, phi] =fmodes (K,M) ; 


wn=sqrt (omega) ; $ Nat Freqs in Rad/sec 
wd=wn* sqrt (1-z.*2);— & Damped Nat Fregs (Rad/sec) 
wd=wd(:,1); 


Qk ke KK RK RK KK KK EKER KK RK KK KEKE KK ER KEK KE 


& Calculate Impulse Response Function 
Dkk KK KR KK KKK RK KKK KK KKK KK KKK KK EKER KK EK 


* 4 


IRF=zeros (1,numstep+1) ; % Initialize IRF vector 
for i=l:numstepil % Set up i-calc-loop 

t= (11) *ae- 

£or LT=Lsdor % Start r-loop for IRF summation 
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if wn(r)==0 % Calculate RBM IRF 


PRE tt) S Re () pn (4) pn.) eG 
else % Calculate Non-RBM IRF 


IRF (1)=IRF(1)+phi(1,r)*phi(1,r) *... 
exp (-z (xr) *wn(r) *t)*sin(wd(r)*t) /wd(r); 


end 
end 6 End r-loop 
| End i-loop 


oe 


end 


ZREKKKKKKKEKEKKKEKKKEKEKKEKKEKEKEKEKKEKKKEKKKKKEKEKKEKK 


%$ Initialize and conduct Synthesis 
KEKKKKKKKKKE KEKE KEK KKKKKE KEK KKK KE KKK KKK KKK KK EK 


oe 


force=ones (numstept+1,1); % Begin with an init f vec 
dif=100; % Initialize difference 
tol=le-2; % Define error tolerance 

t= baat cells % Define Time vector 
y=Y0*sin (W*t) ; % Define Forcing Function 
t0=clock; 


for ict=1:300 


oP 


Begin synth iteration, 
% start ict-loop 

cc=conv (IRF', force); 

x=-dt*ecc(1l:numstep+1,1); 


if 2CES=2 % Begin error determination 
[11,j3}]=max (abs (xlast1l-x)); 
Gif=max (abs (xlast1(jj)-x(3j3))); 

end 


xlasti=x; 


LE - ALE<tol 


oe 


Check error vs tolerance 
% Plot results if acceptable 
disp('Breaking'); | 
disp(sprintf£('Iterations to convergence... 
eo1. 7 Let) ) 3 
figure (1) ,ploc(t;s;"o';t,*lasti,"2r') 


xlabel ('Time(sec)'),ylabel('Displacement (in) ') 
title('Last(red) vs New(green) ') 
break 
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end 
[force] =1500* (x-y’); 


end % End i-loop 
timei=etime (clock, tQ); 


eK KK EK KK KKK KK KK KEKE KK KEK KKK KEK KEKE RK K KKK KKK KEKE KREKEEE 


de 


Perform exact solution using ODE23 
Pe eee reer Cr ener er eer Tee re 


de 


oP? 


Initial displacement/velocity 
vector for ODE23 
Time vector for ODE23 


ae 


xi=[0 0 0 0]; 


de 


oe 


tspan=([tisdt: tli; 
tl=clock; 

[Tode , Xode] =ode23 (’vibe2a’,tspan,x1) ; 
time2=etime (clock, til); 


figure (2) 
plot (t,xlast1,‘r’,Tode, Xode(:,2),'x'} 


mabe 








Ne ee ee 














APPENDIX B: LINEAR AND NONLINEAR SYNTHESIS COMPUTER CODES 


USED FOR EXAMPLES 3-1 and 3-2 


The following is a brief description of the MATLAB computer codes used to 


perform the calculations necessary for Examples 3-1 and 3-2. The code for both 


examples is nearly identical with a minor modification to replace the linear stiffener with 


a nonlinear stiffener in Example 3-2. Both examples are of a large plate-mass system 


(free-free) which is modified by connection to ground at the four plate corner points. 


NLISOTEST.m — reads model data from a NASTRAN PCH data file via 
subroutine (fReadNASModesPCH.m). Applies data from NASTRAN file to 
another subroutine (fIRF.m) which calculates the Impulse Response Functions 
for the points of interest/modification in the model. Also generates the 
forcing function via another subroutine ({BlastForcing.m) and feeds all 
relevant data to the synthesis subroutine (fSynth2.m). 

fSynth2.m — first uses data provided from NLISOTEST.m to aaah the 
post-modification transient responses of the four plate corner (connection) 
points. The transient response of each corner is calculated separately and the 
appropriate IRF is extract from the I RF matrix using a subroutine 
(fSymmetricStore.m). Transient response is obtained through an 
iterative/convergence algorithm which uses an initial guess for the force 
vector to generate an initial displacement vector. This is applied to the actual 
isolation system, which produces the next round of force and displacement 
vectors and convergence calculations are conducted. In Example 3-1 a linear 
stiffener is used, while in Example 3-2 a nonlinear stiffener is used 
(fNonlinearSpring.m). Based on the transient responses of these corners, the 
program calculates the transient response of the installed equipment using the 


equipment IRF data. 
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fIRF.m — uses data from NASTRAN model and defined points of 
interest/connection to develop unique Impulse Response Functions and stores 
them in a symmetric column. Allows for rigid body modes by assuming any 
natural frequency less than le-3 (rad/sec) is a rigid body mode. 
fBlastForcing.m — can be used to generate one of three (sine, step or blast) 
output forcing functions based on input data. A blast function is used in this 
thesis. This subroutine produces displacement vs time histories for all three 
types. Also, for sine and blast functions, velocity vs time histories are 
generated. 

fSymmetricStore.m — is used to extract the appropriate IRF from the IRF 
matrix when calculate the transient response of the four plate corners. 
fNonlinearSpring.m — is used as the force-displacement characteristic of the 
nonlinear stiffeners installed in Example 3-2. This subroutine uses a table- 
lookup method for determining the forces generated by a series of 


displacements based on the data points contained in the program. 
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NLISOTEST.m 


global x 


zeta = 0.02; 


oe 


Se KKK KK KK KKK KKK KKK KK KKK KKK KKKKEKKEKEKKKEKEKE 


Read model data from NASTRAN PCH File 


Qk kK KK RAKE KK RK KKK RK ERK KK KKK KKK KKK KER ERE KE 


de 


[lam, phi,grid_lbls]= fReadNASModesPCH (’P16900Modes.pch’); 


2 


= sqrtilam); SCalculate Nat Fregs 
wn (1:9); sRetain only 1°* 9 nat fregs 


wml 


% Time Step: 


GPa eee 

Start. = 020; % Start Time 

del_t = le-4; % Time step size 

end_t 04; % End time 

time = [start_t:del_t:end_t]’; % Time points 
nstep = length(time) ; & No. Time points 
“b set = [3:6:27]; % Define points of importance 


Qe KKK KERR RK KKK RK KKK KKK KKK RR KK KKK KKK EKER KS 


% Calculate system IRF 


kk KERR KEK KKK KEK KK KEK AKER KEK KKK KKK KK ERS 


[irf] = fIRF (wn, zeta,phi, b_set, time); 


Rea K KKK KEK EK KKK KR KEK KER KEK KK KK RK KEKE RK EKER . 
% 


Generate forcing function (BLAST) 
eK KKK KK KKK RK KK KK KKK KKK KEKE KKK KEKE KEKKEKEKEEEESR 


oe 


Fo = 1; 
type = ‘blst’; % sine 
plotit = 0; 


[y_of_t,ydot] = fBlastForcing(Fo,time, type, HLOCLTE)s 
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QRKKKKKKEEKEKKEKEKK KER KEKEKKEKKEKKEREEKKKRKEKRKEKKEKEKEKKERKEEKK 


% Start synthesis, time synthesis 


DEREK KKEKKEKEKKEKKEKREKKEKKEKKEKEK KER KEKEKEKKEEKKKEKRKERE 


tO=clock % Start time for clock 
iter _limit = 400; Max convergence iterations 


[x,£,icnt_iter]=fSynth2 (time, length(b_set) .... 
Let; yor tere, 11m): 


regqtime=etime(clock,t0) % Time reqd for synthesis 
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fSynth2.m 


function (x, f,cnt..1ter)]=Synth2 (time; num-bset, irt... 


Vv .0r t,1Lter. amie) 


global x iter 


of oP 


oe 


dP dP dP oP 


dP dP oP oP 


oP o& dP dO 


JP oP dP dP AP DP DO cP ADO GP HP oO 


dP dP de dO 


oe 


of co oP 


de 


Usage: [x,2,i¢ent. 1ter)|=ESynth (cime,b Set, ir, y-or tc); 


input Arguments: 


time: Sample points, size(time) = nstep x 1 


num _bset: Number of base excitation coordinates. 
The size of the synthesis problem is 
length(time) * num_bset 


pi apie Matrix containing the impulse response 
functions stored in symmetric column storage 
(see fSymmetricStore.m). 
The number of columns=num_bset’*... 
(num_bset+1)/2 


y_of_t: Base displacement time history. Applied to 
all b_set points. 


iter _limit: Limit on number of iterations. Solution 


is returned with a warning. 


Output Arguments: 


at gh Menge gt gt gt Sgt gt Sag Pt Cagt Pia Ting Pint Pig Pet Mngp 


x: Synthesized response of b_set coordinates. 
Size(x) = nstep x num_bset 

£: Synthesized forces acting at b_set 
coordinates. 
size(f) = nstep x num_bset 
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de 


oe 


Coordinate Sets: 


Pt a et Fe at i et Png Map Pine gt gt Fat Cnt as Tal 


ae dO 


oe 


* Coordinate labels refer to FE model coordinates 
Row vectors which are lists of the 
coordinates in the four sets: 


de dO ce 


fo) 


% where {b_ set} = vector of base displacement 
% coordinate 
% {m_set} = vector of modification 
% | coordinate 
{c_set} = vector of coupling coordinate 


de dO 


{i set} vector of internal coordinate 


fa) 
1] 


oe 


Currently, the labels are not used in this 





% function. 

2. 

(*) 

% 

nstep = length(time) ; % No of steps 

del_t = time(2) - time(1); % Time step size 

x last = ones(nstep,num_bset); % Initialize x_last 

f£ = ones (nstep,num_bset); % Initialize force 
converge_tolerance = 1.0; % Convergence criteria 
converge_check = 100; % Init Convergence Value 
1ent.1ter = 0; % Init Iteration Counter 
Dok kk RR RR KK RK RR KEK KERR KK KK KKK RK KKK RK ERE KEKE KK 
% Synthesis Iteration for transient response of plate 
6 corners 

Qk KEKE KKK KK KAKA E KEK HEH RHEE EK EEK KEKE KEE KK KK KKK K EK 


while converge_check > converge_tolerance 


x = zeros(nstep,num_bset) ; 
xdot = zeros(nstep,num_bset) ; 
Tent Ater = acnt-i1ter + 1; 


for icnt_bset_rows = 1 : num_bset-l,; 


for icnt_bset_cols = 1 : num_bset,; 
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[sym_col] =fSymmetricStore(icnt_bset_rows.... 


icnt_bset_cols,num_bset) ; 


fi COnv: + =—Cony (17 FCs, Sym Cot) 72% 
£(:,icnt_bset_cols)); 


x(:,icnt_bset_rows) = x(:,icnt_bset_rows)... 


- del_t * h_conv_f(li:nstep,1); 


end 
ydot=fdot (y_of_t,del_t)'; 


xdot(:,icnt_bset_rows) =fdot... 
(x(:,icent_bset_rows) ,del_t)'; 


£(:,icnt_bset_rows) =fNonlinearSpring... 
(x(:,icnt_bset_rows) - VT On. 0 hs 
.02* (xdot (:,icnt_bset_rows) -ydot) ; 


end 


Qk KK KR KK KKK RK KKK KKK KER KKK KKK RK KKK KEK KE 


& Conduct convergence check and verify 


% Max iterations have not been reached 
Cea KKK KK KKK KKK KEK KKK KR KEK KK RK RRA RE REE KE 


converge_check = fConvergeCheck(x,x_last) ; 
x last = x; 
if icnt_iter > iter_limit; 
fExceed_Iter(icnt_iter) 
break; 
end 


end 


£in=(5 9 12.14]; 


67 





KEKKKKEKKEKEKEKEKKKE KE KKEKEKKEKEKRKKEKEKKKEKKEKKKEKKEKEEESE 


Calculate Transient Response of Equipment 
KK KKK KKK KEK KKK KKK KR KEKE KKK KKK KEK KK KEKE KKK REE 


of oO 


fe\®) 


for 1-1 + lengeht{ fin) 
h- conv. fi <22 )=conv (irk (stints) ) ek 4,1.) )3 
x(:,5)=x(:,5)-del_t*h_conv_f(1l:nstep,1); 
end 
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£LIRF.m 


fIRF(wn, zeta, phi, dof, t) 


function -LLRF “cym.col | 


% Usage: [IRF_sym_col] = fIRF(wn, zeta, phi, dof, t); 


® Function uses all modes passed in. 

® Function will generate all unique input-output pairs 
% defined by dof. 

& Function will store them in symmetric column 

% storage. 

% See f£SymmetricStore.m) 


S IRF are evaluated at the time points in the vector 
% ts 

% wn is in rad/sec 

% mass normalized phi assumed 

& zeta is a scalar, and is applied to all modes. Zero 
$ 1s OK. 


% If wn < le-3, rigid body mode is assumed. 


num_modes=length (wn) ; 

num _dof=length (dof) ; 

num_sym_col=num_dof*num_dof+1)/2; % Number of columns 
IRF_sym_col=zeros(length(t),num_sym_col) ; 

S$Initialize matrix 

mode_irf=zeros(length(t),1); 

icnt_sym_cols=0; 


% Will end up being num_sym_col 


for Lent row der = 1. = numdot; 


for icnt_col_dof = icnt_row_dof : num_dof; 
icnt_sym_cols = icnt_sym_cols + 1; 
for icnt_modes = 1 : num_modes; 


if abs(wn(icnt_modes)) > 0.1; 


Then elastic mode 
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wd = wn(icnt_modes) * sqrt(l1 - zeta%2); 
mode_irf=exp (-zeta*wn(icnt_modes) *t)... 
.*sin(wd * t) / wd; 
elseif abs(wn(icnt_modes)) <= 0.1; 
SRigid body mode 
mode_irf = t; 


end; % End if abs(wn) 


row_dof = dof (icnt_row_dof); 
col_dof dof (icnt_col_dof); 


IRF_sym_col(:,icent_sym_cols)=IRF_sym_col (:,.. 
icnt_sym_cols)+phi (row_dof, icnt_modes) *... 
phi (col_dof, icnt_modes) *mode_irf; 


end; % End icnt_modes 
end; % End icnt_col_dof 
end; | % End icnt_row_dof 


% End function fIRF.m 
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FblastForcing.m 


function[f_of _t, fdot]=fBlastForcing (Fo,time,type,plotit) ; 


Usage: [f_of_t, fdot]=fBlastForcing 
% : (Fo, time,type,plotit) ; 


&% Choices: sine blst step 


% type='step’ STRING Variable 


If use ‘sine’, fdot also returned. 


% This function returns a forcing function 
® which is a "blast" function 


% F(t) =Fo* (exp (-at) -exp(-bt) ) 


% Where a and b are constants which shape the blast, 
and Fo is the amplitude of the blast 


2 The variable "plotit" is a switch which if set=1 
® will cause the f(t) to be plotted, if set to 
% anything else will not plot. 


af type==' bist’; 
Sdisp(’Blast forcing used...’);_ 
a=100.0; 
b=300.0; 
f of t=Fo* (exp (-a*time) -exp(-b*time) ); 
fdot=Fo* (-a*exp (-a*time) +b*exp (-b*time) ) ; 


elseif type==’step’ 
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S$disp(’Step forcing used...'); 
f_of t=Fo*ones(size(time) ); 


fdot=[]; 

elseif type==‘sine’ 
$disp(’Sine forcing used...’); 
W=5; %Hz 


f_ of _t=Fo*sin(2*pi*w*time) ; 
fdot=Fo* (2*pi*W) *cos (2*pi*w*time) ; 
end 


LE pPlotit==1;> 
figure (gcf+1) 
if type=='sine’ ; | 
DLlot(time, t..0of t;time; fdot)+ grid 
else 
plot (time, f_of_t);grid 
end 
pause (2) 
end 
6 END FUNCTION 


[o) 
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fSymmetricStore.m 


function [index] =fSymmetricStore(row,col,n) ; 


dP dP ch dO BO BO BO oP BHO HH CO 


oP 


dP dO dP oP 


de 


oe 


d& db oP 


dP dP 


de dP dO dP ce oO 


oe 


Usage: [index]=fSymmetricStore(symcol,n) ; 
Or: [index]=fSymmetricStore(row,col,n); 


This function returns the row and column indices 
corresponding to the column number of a matrix of 
size n where only its upper triangle is stored 
(i.e. col>=row); 


~—-—-Or--- 


The function returns the column number corresponding 
to upper triangle symmetric storage, given a row/col 
index pair. 


Example: A 4x4 matrix is shown. Element (2,3) has 
a "symcol" of 6 


4 2 3 4| | 

| 5 6 7| stored columnwise as [123456789 
| g 9| 10] 

| 10| 


Written by J.H. Gordis, Ph.D. 
Rev 1: 9/15/98 


Input is symcol and n -> Output is (row,col): 


Pe ce Pee ad at Rint eat iat at Pd Ph int at Prat Prat Oat Prat Pind Prat ad Mat Pat Pat Pat at ad Peal Pact Pap Pt Prat Pet nat Pend Paap Rat Peat Mat Pap Peat up Meet Pret Peet Peat at Cat et Naf Oral Cot 


if nargin == 2 & row <=col*(col+1)/2; 


rowout=ceil (row/col) ; 
colout=(rowout-1)+rem(rowout,col) ; 
index=[f[rowout colout] ; 


elseif nargin ==2 & symcol >n*(n+1)/2 


disp("’) 

disp(‘Error in fSymmetricStore’ ) 

disp('Symmetric storage index exceeds max value... 
for matrix size "n". ') 
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disp (’*) 
elseif nargin==3&row<=n&col<=n&col>=row 


SInput is (row,col,n) -> Output is symmetric column 
% number: 


6 Number of elements contained in j rows of upper 

% triangle: 

% 3*n-sum(1:j-1) 

% Add up all elements in rows 1:row-1 (j=row-1). Then 
& add additional elements for col. 


index= (row-1) *n-sum(1:row-2)+(col-row+1) ; 
elseif nargin == 3&row<=né&col<=n&col<row 


temp=row; 

row=col; 

col=temp; 

index= (row-1) *n-sum(1:row-2)+(col-row+1) ; 


de 


disp(’’) 

disp('Error in fSymmetricStore.’) 

disp{'Row>Col: Indices refer to element in lower... 
triangle.') | 

disp(' ') 

elseif nargin==3&row>n|col>n&col>=row 

disp('') 

disp('Error in fSymmetricStore') 


de 


de dP oO 


disp ('Row and/or Col indices exceed matrix size "n". 


disp('') 
end 
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fNonlinearSpring.m 


function [fofx]=fNonlinearSpring (xin, klinfit,plotit) ; 


KK KKK KKK KKK KKK RK KKK KKK KKK KKK KEKKKEKKEKREEE 


de 


& Define Spring Force VS Deflection 
KKK KK KKK KR KKK KKK RR KKK KKK KERR KKK KE EEKREEKE 


oe 


xvst= [0.0 0:0; 
20020; 


™ 


0 
400.0 
600.0; 
780.0 
960.0 
1140. 
1320. 
1500. 
1600. 
16502 
1700. 
1740. 
1750. 
1760. 
L765 
1770. 

. 1770: 
10 0: 1770.0; 
100.0 277020; 
1000.0 1780.0; 
1900.0 1790.0; 
2800.0 1800.0; 
3700.0 1810. 
4600.0 1820. 

8) 
0 
0 
0 
0 


PERE PBEPEPERPOODOOOCCCOO 
ONUPBPWNHEFOWDOIHDUBRWNHE 


OD 0.0 0:-O oOo OO 0 0 © 'O Oo oO 


=e 


“=e 


5500.0 Les0s 
6400.0 1840. 
TI00<'0' 650% 
, 8200.0 1860. 
~°9100.0 1870. 
10000.01880. 


=e 


~e e = 


MDrOdDOD OOO CO 


L—d 


e 
f 


num_pts=length(xvsf) ; 


i ie 














ie) 


% Reflect points across x-axis 
[x,£]=£XYreflect (xvsf(:,1),xvsf(:,2)); 


ce Ki He te SO; 
FIANEVCHkliniite*xvst( 1): 
[xX ,£lintit |/=txYreéfiect (xvst (=,1);flintit)-; 
FOLf=f—-fLiniitc. 
else 
fdif=zeros(size(f))’; 
flinfit=zeros(size(f))’; 
end 


£OLx=1INnterp L(x, foci): % Interpolate force 
% via table lookup 


BREKKKEKEEKKEERKEKRKEKREKKEEKEKREKKEKEKEKKK 


Plotting functions 
ZREEKKKKKEERKEKKEKEKKEKKEKEREKEREEKEKEKKEEKE 
if nargin==3; 
plot(x,f,'r’);grid;figure(gcf) 
title(’Nonlinear Spring Force VS Displacement’ ) 
xlabel(’Deflection (in) ’) 
ylabel (‘Force (1lbf£)’) 
end 
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APPENDIX C: OPTIMIZATION COMPUTER CODES 


The following is a brief description of the MATLAB computer codes used to 


perform the optimization routine. Many of the subroutines used in the optimization 


routine are identical to those presented in Appendix B, and are not presented here. Code 


presented here solves the nonlinear stiffener modification only, but may be easily 


modified to solve a linear stiffener modification problem. 


fMainOpt.m — performs all of the functions of NLISOTEST.m presented in 
Appendix B. Also defines the starting point for the optimization routine and 
determines the time required for completion. This program uses CONSTR.m 
and calls Optfun.m to perform the optimization. 

Optfun.m — invokes the synthesis subroutine (Synth2opt.m) to generate model 
transient response data. The program then uses the transient response data to 
calculate values which are needed for objective function and constraint 
function evaluation. 

Synth2opt.m - similar to the subroutine fSynth2.m presented in Appendix B. 
The primary difference is that this function generates a nonlinear spring force- 
displacment characteristic base on design variables. Since a table-lookup 
scheme (such as the one used in fNonlinearSpring.m) cannot be used for this 
problem, a decision algorithm is used to determine the forces generated by the 
spring base on the displacement time history. First, the transient responses of 
the four plate corners are synthesized. Second, the transient response of the 
equipment is calculated using plate corner response data and the appropriate 


IRF’s. 
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fMainOpt.m 


KKK KKK KEK KEK KKK KKEKKKKAEKEEEE 


MAIN OPTIMIZATION PROGRAM 


KKK KKK KEK KEK KK KKK KEKKKEKEKSE 


dP oO 


oe 


clear all; 
close all; 
Cilc* 


global irf y_of_t xout iter kkmod xx ydot force 
global forces ff maxac rlbs rlcmp energ 
Zeta=0.03: : 


Se KK RR KK KK KK KK RK KKK KKK KEKE KR KK KKK KKK KKK KEKE KKERER 


Read NASTRAN data file 


eK KK KK KEKE KKK KR RE RE KKK EEE KEKE KKK KKK KEKE KK KEKE KK KK 


de oO 


de 


(lam, phi, grid_lbls]=fReadNASModesPCH (’P16900Mo.pch’); 


wn=sqrt (lam) ; : % Calc Natural Fregs (Rad/sec) 
wn=wn (1:9); % Retain ist 9 Nat Freqs 


ge KK KK KK KK KE KK KE KK KER KEKE KKK KKK KKK KKK KEK KKK KE 


de 


ree 


Generate time vector and define # of steps 
eK KK RK RK RK RK RK KKK KKK KKK EK KEKE KEKE KER EK KEKE KE KK KKK 


fone 


start_t=0.0; 
del t=l1e-4; 
end t=0.04; 


time=[{start_t:del_t:end_t]’; 
nstep=length (time) ; 


eK KK KR KE RE KR KK KR RK KKK KEKE RE KK KEK KEKE KEKE KEKE KEKE KEKE 


de 


dP 


Define points of interest and connection points 
Calculate IRF for each 


2 RR ak RR RR Ra RR RR RK RK KK RK KR RK KKK 


oe 


oe 


b: set=(326-27); 


[irf]=fIRF (wn, zeta,phi,b_set, time) ; 
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SEERA AREA RAR HERE AR CK RK RS RRR AR eR 


% Define forcing function 
KKK KK KK KEK KKK KKK RK KKK KKK KKK KKK EKER KEKE KEKEKEKEEKRE 


de 


Fo=1; 
type="bist’; 
plotit=0; 


[y_of_t, ydot]=fBlastForcing(Fo,time,type,plotit) ; 


iter_limit=400; % Synthesis iteration limit 
1Ler=0: % Optimization counter 
LEEK ELLER ELEL SESS SAK AERA RA EARA SAAR RRS 


ae 


Define starting point and upper/lower bounds 


for optimization routine design variables 
KK KK KK KK RK KKK KEKE KEKE KEKE KKK RRR KKK KEKE KKKEKKEKEKRER 


de 


* 


fee 


kmod0=[1.00 9.65 .06 .0705 .5]; 
kmodub=[1..5 -10.0 .07 .071 5.013 
kmodib=[.5 7.21 .045 .0701 0.0]; 


KK KKK KKK KEK KKK KK KEK K KKK KKK KK KK KEKE KRE KK KK KKKEKEK 


% Define Options to be used in CONSTR.m 


Be KK KKK KK KK KK KKK KK KEKE KKK KKK KKK KKK KE KKK KKK KEKE 


oo 


de 


options (1)=1; 
options (14) =1500; 


Fem KKK KE KK KKK KKK KK KEE KE KK KEKE KEK KKK RE KKKEKREKKEEERER 


oe 


de 


Begin and determine time required to perform 
Synthesis. : 


KKK KKK KK KK RK KK KK KEKE KR KRE KK KEKE KEKKKEKRKEKKEKKKEKEKEKEKERE 


de 


de 


t0O=clock; 
[kmod]=constr('Optfun’ ,kmod0, options, kmodlb, kmodub) 
reqtime=etime (clock, t0) 


save RUNnew4damp2 % Save data for test 
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function[£,g]=Optfun (kmod) 


global irf y_of_t xout iter kkmod xx ydot force ff 
global maxac rlbs rlcmp energ 


KK RK KEK KK KE KEKE KKK KEK KEK KKK KKK KKK KEKE KEKE KREKRKKKKAKKEKSE 


relied 


oe 


Generate time vector and define # of steps 
KEKE KK KEKE KKK KE KK KKK KEKE KKK KKK KKK KK KE KKK KEKEKKREKE KEKE 


do 


start_t=0.0; 
del _t=le-4; 
end _t=0.04; 


timne=[starttsdel.t send. tt]; 
nstep=length (time) ; 


b_set=[3:6:27]; 
iter _limit=400; % Synthesis iteration limit 
Teer=1ctert+i; % Optimization iteration counter 


KREKKKKEKK KE KK KEKE KK KEKE KEKE KEKKEKKKEKEKEKEKAEKEKEE 


oe 


de 


Call synthesis routine 
RkKKK KK KEKK KK RE KKK KKK KEKE KEK KR KKK KK RK KK KK 


[xout, force, kmod]=Synth2opt (time, length (b_set) .... 
TEV 208. 0.2 Cer 11m emod):; 


KK KEK KEKE KKK KKK KEKE KKK KEKE KEK KEK KEKE KK KEKE KE KEKE KKKEKEKEKERKSE 


de dP 


Calculation of values to be used in constraints 


and relevant optimization tracking data 
KKK KKK KE KK KKK KKK KK KEK KEKE KKK KR KEK KKK KKK KEKEK KEK KKK 


de 


ae 


KEKE KK RK KKK KEK RK KK KEKE KEKE KK KEK KKK KEKE KK KEKEEKEKKEEKESE 


oe 


oP 


xddot--equipment acceleration; acc--max equip accel 


in g; maxac--max equip accel experienced for each 


Contigquratwon 
KK KK KR KK KK KK KKK KKK KEK KEKE KKK KKK KEKE KKK KEKE KKEKKKEKEKKKAEEESE 


oP dO 


ae 


xddot=fddot (xout(:,5),del_t); 
acc=max (abs (xddot) )/386.4; 
maxac (iter) =max (abs (xddot) )/386.4; 
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KK KKK KK KKK KKK KKK KEK KEKE KKK KKK KKK KKK KKK KK KEKE KEKEKKEKEEKER 


kkmod--track design variables for each configuration 
xx--track equip absolute displacement " " " 
relbase--track base relative displacement 
relcemp--track equip relative displacement 
s(1,2)--track stiffener characteristics " '" 
energ--track equip displacement away from 


equilibrium " " 
KK KKK KKK KKK KKK KKK KKK KKK KEK KEK KKK KEKE K KERR KEKE KRKKEKESE 


de od? 


oe 


bi it 


oe 


W ci bd 


dP dP dO dP 


oe 


kkmod (iter, :)=kmod;. 

xx (2, 2 ter) =xout(s:,5)4 
relbase(:,iter)=(xout(:,1)-y_of_t); 
relcomp(:,iter)=(xout(:,5)-xout(:,1)); 
rlbs=max (abs (relbase(:,iter))); 
ricmp=max (abs (relcomp(:,iter))); 

s(1)=kmod(3) /kmod(1)*1e4; 
s(2)=(kmod(4)-kmod (3) )/ (kmod(2)-kmod(1)) *le4; 
energ (iter) =0.5*trapz(xout(:,5).%2); 


KK KK KK KK KKK KE KKK KE KKK KEK KEK KKK KEK KR KKK KEK KE KKKEKKKEKKEKEEKEK 


fohes 


Define objective function as minimization of max 


accel of equipment 
eK KK KK KKK KK KK KKK KEKE KKK KKK KK KKK KKK KE KEKE KEKE EKKEKKEKSE 


de dP 


de 


f=maxac (iter) ; 
£f (1 Ler) =r 
xmax=max (abs (xout(:,5))); 


KKK KEKE KKK KKKKEKEKRKKEEKESK 


oe 


Define constraints 
KKK KKK KEK KK KEKKEKEKEKEKKEKEKK 


oe 


de 


g(1)=rlbs-.50; 
g(2)=rlemp=. 50; 
g(3)=xmax-1.0; 
g(4)=acc-16; 


% End bdfun2acom.m 


o 
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Synth2opt.m 


function[xout, force, kmod]=Synth2opt (time,num_bset,irf,... 
y_of_t,iter_limit,kmod) ; 


global xout ydot forces 


KEKKKKEKRKEEKRKER ER KEKKEK KEK KER EK KEK KEKE KEK KEEKKE KEKE KEK KEK KEKE KKK KREKEKEEKEK 


de 


de 


Calculate slopes of force-displacement curves based 


on the design variables 
KKK KEKKEKEK KEK EKK KEKE KKK KE KEKE KEKE KKK KEKE KKK KEKE KKK KEKE KKK EES 


oe 


oe 


s (1) =kmod (3) /kmod(1) *1e4; 
s (2) =(kmod (4) -kmod (3) ) / (kmod (2) -kmod (1) ) *1e4; 
kmod (5) =1le-2*kmod (5) ; 


KKK KKK KK KEKE KKK KEKE KEKE KEKE KK KKK KKK KEKE KKK KKK KKK KKK KKEKKEEKSE 


de dO 


Define number of steps (nstep) , time step size (del_t), 
Initialize x_last for convergence calculations and 


Make initial guess for starting force vector (>0) 
KREEKKKEKKK KK KK KK KEK KK KEKE KKK KEK KEKE KKK KK KK KEKE KEKE KEKE KKKEEKEKEESE 


ae oP 


oo 


nstep=length(time) ; 
del_t=time(2)-time(1); 
x_last=ones (nstep,num_bset); 
force= le-9*ones(nstep,num_bset) ; 


converge_tolerance=1; % Initialize convergence requirements 
converge_check=100; 
icnt_iter=0; SInitialize synthesis iter counter 


while converge_check>converge_tolerance 


xout=zeros (nstep,num_bset) ; 
xdot=zeros (nstep,num_bset); 
icnt_iter=icnt_iter+1; 


for icnt_bset_rows=1:num_bset-l1; 


for icnt_bset_cols=1:num_bset-1; 
[symcol]=fSymmetricStore(icnt_bset_rows.,... 

icnt_bset_cols,num_bset) ; 

h_conv_f=conv (irf(:,symcol),... 
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force(:,icnt_bset_cols) ); 
xout(:,icnt _bset_rows) =xout(:,icnt_bset_rows)... 
-del_t*h_conv_f(1l:nstep,1); 


end % End icht_bset_cols 


Skt RRA ARI OR Rae eK IK eR A I RK eR RICKS RR RARE 


% Force calculations for plate 
TALKER LAK ERAAAKA ALA RARER AAAS A RAE AA AREA AAS 


ydot=fdot (y_of_t,del_t)'; 
xdot(:,icnt_bset_rows) =fdot (xout(:,icnt_bset_rows) ,... 
del _t)'; 


KKK KKK KKK KEK KKK KKK KEKE KEKE KEKE KKK KKEKKEKKEKKREKREEEKEK 


Separate response vector into regions of 
force displacment characteristic curves 

for real force calculations. First conduct 
a logic check to verify if response is 
within some bounds. Multiply logic vector 
by original vector to generate separated 
response vectors. Use separated response 
vectors to generate separated force vectors 
and combine these force vectors to obtain 


real/overall force vector. 
fe KKK KK KK KKK IKKE KKK KEK KEKE KEK KKK KKK KKK KEKE KEKE 


dP dh BP GD oH GP BO BH DO CO AO 


dP 


if icnt_bset:rows < 5 


qm(:,icnt_bset_rows) =abs (xout (:,icnt_bset_rows)... 
| -y_of_t)<=kmod (1) ; 
qu(:,icnt_bset_rows) =(xout(:,icnt_bset_rows)... 
=v OF tC) >kmod (1); 
ql(:,icnt_bset_rows) =(xout(:,icnt_bset_rows)... 
~y_of_t)<(-kmod(1)); 
rmm(:,icnt_bset_rows)=gm(:,icnt_bset_rows) .*... 
(xout(:,icnt_bset_rows) -y_of_t) ; 
ru(:,icnt_bset_rows) =qu(:,icnt_bset_rows) .*... 
(xout(:,icnt_bset_rows)-y_of_t) ; 
r1(:,icnt_bset_rows)=ql(:,icnt_bset_rows) .*... 
(xout(:,icnt_bset_rows) -y_of_t); 
fm(:,icnt_bset_rows)=s(1)*(rm(:,icnt_bset_rows) ) ; 


fu(:,icnt_bset_rows)=(1le4*kmod(1)*kmod(3).*qu(:,icnt_bset_rows))... 
+(ru(:,icnt_bset_rows)-kmod(1).*... 
qu(:,icnt_bset_rows)*s(2)); 
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fl(:,icnt_bset_rows) =- (1e4*kmod(1)*kmod(3).*.. 


gl(:,icnt_bset_rows) )+ (rl... 
(:,1lcnt_bset_rows)+kmod(1).*... 
gql(:,icent_bset_rows)*s(2)); 


force(:,icnt_bset_rows) =fm(:,icnt_bset_rows) +... 
fu(:,icnt_bset_rows) +... 
fl(:,icent_bset_rows) +kmod(5) *... 
(xdot(:,1lent_bset_rows) -ydot) ; 
end 6 End if loop (plate force synthesis) 


end % End icnt_bset_rows 


KERRIER RRR RE CRERREEK KEKE K KK RE KEK 


de 


oe 


Check for convergence 
KKK KKK KKK KKK KE KEKE KKK KKK KEKE KKK KEKE KEKE KEKE 


de 


converge_check=fConvergeCheck (xout, x_last); 


x last=xout? 


1f icnt_iter > iter limit; 
fExceed_Iter (icnt_iter) 
break; | 

end 


end % End while loop 


dP dD ce cO 


de 


KKK KEKK KERR KEKEK KEKE KEKE KKK KKE KKK KEKE KEKKKKKKRI EH 


Calculate response of equipment based on 
synthesized response of plate corners and 


equipment impulse response functions | 
KEKE KE KEK KK KEKE KEKE KKK KEKE KEKE KKK KEKE KKK KE KKK KEKEKS 


Fin=[5 9 12 14]; 


For 2H=l<lengtht{ iin) 


hh conv. f (+, 1)=conv (1rt (:, fini) ) , force (<,1))> 
xout (:,5)=xout(:,5)-del_t*h_conv_f(1:nstep,1) ; 


end 
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